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We present a first-principles based multiscale modeling approach to heterogeneous catalysis that 
integrates first-principles kinetic Monte Carlo simulations of the surface reaction chemistry into a 
fluid dynamical treatment of the macro-scale flow structures in the reactor. The approach is applied 
to a stagnation flow field in front of a single-crystal model catalyst, using the CO oxidation at 
RuO2(110) as representative example. Our simulations show how heat and mass transfer effects can 
readily mask the intrinsic reactivity at gas-phase conditions typical for modern in situ experiments. 
For a range of gas-phase conditions we furthermore obtain multiple steady-states that arise solely 
from the coupling of gas-phase transport and surface kinetics. This additional complexity needs 
to be accounted for when aiming to use dedicated in situ experiments to establish an atomic- 
scale understanding of the function of heterogeneous catalysts at technologically relevant gas-phase 
conditions. 
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I. INTRODUCTION 

First-principles kinetic Monte Carlo (lp-kMC) simula- 
tions have evolved into an important tool in the mod- 
eling of heterogeneous catalytic processes 1 . The suc- 
cess of the approach relies on the accurate treatment 
of two central aspects for the reactive surface chem- 
istry: A first-principles description of the involved ele- 
mentary processes and an evaluation of their statistical 
interplay that fully accounts for the correlations, fluctu- 
ations and spatial distributions of the chemicals at the 
catalyst surface. Particularly if suitably combined with 
sensitivity analyses 2 , lp-kMC simulations thus offer the 
prospect of an error-controlled and quantitative microki- 
netic modeling of the surface catalytic function. In par- 
ticular for technologically relevant environments, i.e. at 
near-ambient reactant pressures and elevated tempera- 
tures with concomitant higher product formation rates, 
a third aspect comes into play that needs to be accounted 
for to properly describe the observable catalytic conver- 
sions. This is the flow of heat and mass in the given reac- 
tor geometry, for instance the transport of formed prod- 
ucts away from the active surface and how efficiently the 
large amount of heat generated by the exothermic surface 
reactions can dissipate into the system. 

Corresponding macro-scale flow structures are suitably 
described at the continuum-level by the transient Navier- 
Stokes equations together with energy and species gov- 
erning equations. The methodological objective of the 
present work is then to couple lp-kMC into such a fluid 
dynamical (FD) framework, thereby augmenting the ac- 
curate treatment of the reactive surface chemistry pro- 
vided by the prior technique with the capability to ac- 
count for heat and mass transport effects. With a brief 
account of the main results already given in ref. 3, we 
focus here in particular on a detailed description of this 
methodology. While the presented approach can read- 



ily be coupled with any computational fluid dynamics 
(CFD) software enabling the treatment of arbitrary reac- 
tor geometries, we develop it in the following specifically 
for a stagnation flow field in front of a flat-faced model 
catalyst. We argue that this is a suitable, though ad- 
mittedly simplified reactor geometry to qualify transport 
effects in modern in situ studies aiming at an atomic-scale 
understanding of the catalytic function of single crystals 
in technologically relevant environments 4 . 

The focus of such studies lies on possible differences 
in the surface chemistry compared to operation in ultra- 
high vacuum (UHV), where the function of model cata- 
lysts has been extensively studied in the past. In order 
to discern corresponding so-called "pressure gap" effects, 
it is important to assess if heat and mass transfer effects 
noticeably mask the true intrinsic reactivity in the in situ 
conditions. Particularly for the often studied CO oxida- 
tion reaction at late transition metal catalysts there are 
good reasons to suspect that transport limitations might 
not be negligible. Typically reported activities indicate a 
high rate of mass conversion at the surface concomitant 
with a large heat release. Using the established lp-kMC 
model for the CO oxidation at RuO2(110) 5 ' 6 as a rep- 
resentative example we illustrate with the coupled lp- 
kMC-FD approach that the peculiarities of the single- 
crystal reactor geometry lead indeed readily to heat dis- 
sipation and mass transport limitations that severely af- 
fect the observable catalytic function. Key factors are 
the degree of heat conduction at the backside of the thin 
single-crystal, and the propensity to build-up a product 
boundary layer above the flat-faced surface. Obivously, 
such reactor-dependent effects need to be disentangled, 
understood and controlled when aiming to compare data 
obtained by different experimental setups, and when aim- 
ing to conclude on the actual surface chemistry at tech- 
nologically relevant gas-phase conditions. 
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FIG. 1: Schematic view of the stagnation flow geometry: The 
gas streams from the inlet (dashed line) towards the flat-faced 
model catalyst of thickness d and positioned at a distance L 
away from the inlet. At the inlet the partial pressures pjf 1 , 
temperature T ml and axial velocity u ml are controlled; the 
radial velocity v lnl is zero. At the backside of the catalyst 
either the temperature J lback is controlled (leading here to 
the isothermal limit) or the heat flux <? S oiid,z is suppressed 
(leading to the adiabatic limit). 



II. MACRO-SCALE FLOW STRUCTURES: 
CONTINUUM FLUID DYNAMICS 

We develop our approach for a simple reactor geome- 
try suitable to discuss heat and mass transfer effects at 
a flat-faced model catalyst. In the so-called stagnation 
flow geometry 7 shown in Fig. I the gas mixture enters 
through an inlet at a macroscopic distance L away from 
the surface. At this inlet the flow is directed towards the 
catalyst surface, with no variation of gas composition, 
velocity, and temperature in the direction perpendicular 
to the flow. As schematically indicated in Fig. 1 the 
advantage of such a geometry is that it results in an ax- 
isymmetric flow profile. Neglecting edge effects, i.e. the 
finite lateral extension of the model catalyst, this flow 
profile can effectively be described by a one-dimensional 
boundary-value problem. As further detailed below, this 
eases the analysis of the influence of the reactor setup 
significantly and allows to extract the relevant physics 
without being riddled by algebra and numerics. 

However, it is not only this symmetry imposed simpli- 
fication which makes the stagnation flow setup appealing 
or better its realization desirable. In the spirit of the 



Surface Science approach to heterogeneous catalysis, the 
fundamental objective of in situ studies of single crys- 
tals is to obtain insight into their intrinsic reactivity in 
near-ambient environments, for example as function of 
temperature T and reactant partial pressures p a . For 
this to be well-defined, a central prerequisite is that all, 
or at least a dominant fraction of the active surface sees 
the same gas phase. If one considers e.g. a flow geometry 
where the stream of reactants would approach the surface 
from the side - thus in some sense an opposite scenario to 
the here discussed perpendicular stagnation flow - then 
this is clearly not the case. Due to the on-going conver- 
sion of reactants into products the gas-phase composi- 
tion close to the surface would gradually change across 
the lateral extension of the single crystal. If there are 
non-negligible heat transfer limitations, this goes hand 
in hand with a non-uniform temperature profile paral- 
lel to the surface. Under such conditions, making a de- 
fined assignment of observed turnovers to specific pres- 
sures and temperatures becomes essentially intractable. 
In contrast, in the stagnation flow geometry at least the 
entire center part of the active surface sees the same gas 
phase, thereby facilitating the analysis. 

In the stagnation flow setup shown in Fig. 1 the most 
relevant spatial coordinate is thus the direction z per- 
pendicular to the catalyst surface. In the axisymmetric 
problem we denote the other, radial coordinate with r. 
Using the inlet height as zero reference for z, the catalyst 
surface is then at z = L and with a thickness d of the 
single crystal its backside is located at z = L + d. In 
a continuum mechanical description the system is corre- 
spondingly characterized by two spatial regions, the flow 
chamber (0 < z < L) and the sample (L < z < L + d), as 
well as three important interfaces, the inlet (z = 0), the 
surface (z = L) and the catalyst backside {z = L + d). 

In the following subsections we will successively de- 
scribe the modeling carried out for each of these re- 
gions and interfaces. In this manuscript this modeling 
aims at a description of steady-state operation, in which 
chemicals get converted at the active surface at a sta- 
ble, time-independent rate. This rate is the turnover fre- 
quency (TOF) in units of molecules per surface area and 
time. We note that this time-independent formulation is 
a convenience, not a necessity. In fact, in particular the 
coupling scheme integrating the surface chemistry into 
the FD environment is also applicable to transient prob- 
lems and we will briefly mention routes in this direction 
throughout the manuscript. In the same spirit it is clear 
that the approach is, of course, not restricted to the sim- 
ple CO oxidation reaction, on which we will concentrate 
from now on for clarity. 



A. Interface I: Inlet 

At the inlet the gas flow is fully controlled. For the 
FD description this defines boundary values for the tem- 
perature T(z = 0) = T , the partial pressures p a (z = 
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0) = pjf (with a =0 2 , CO, C0 2 ), the total pressure 
pini _ ^p™^ and the axial velocity u(z = 0) = u . 
We set the radial velocity to v (z = 0) = as is the 
case in stagnation flow reactors with a sieve-like show- 
erhead as inlet 7 , and we restrict our attention to flow 
situations with no circumferential motion. Instead of 
the partial pressures it is more common to use the to- 
tal density p and mass fractions Y a as independent fields. 
For the present catalytic context the conversion is readily 
achieved through the ideal gas law 



-pk B T 



(1) 



where m a is the mass of species a and fee is the Boltz- 
mann constant. Similarly, we will see below that it is 
more useful to consider the scaled radial velocity V = v/r 
in the stagnation flow equations. 



Species mass 



Axial momentum 



du 



pu- 



dYa 
dz 



dp + 4 d 
^ dz dz 3 dz 



dz 



du ' 
p— - pV 
dz 



n dV 
2p— 

dz 



Radial momentum 
dV 



pu 



dz 



PV 2 



-A r 



dV 

dz V ^ dz 



Internal energy 



dT 



pc p i 



d dT 

dz dz 



dT 

^c p<aJa<z — 



(3) 



(4) 



(5) 



(6) 



B. Region I: Flow chamber 

1. Flow equations 

The continuum mechanical description of the heat and 
mass transport in the flow chamber is based on balance 
equations for total mass, species mass, momentum and 
internal energy. For the present purposes, the general 
formulation is much simplified by treating the gas phase 
as a Newtonian fluid with vanishing bulk viscosity, by 
the absence of relevant gravitational forces and by the 
absence of significant gas-phase chemical conversions in 
the context of low temperature CO oxidation (i.e. the 
reactions are confined to the catalyst surface). Another 
significant simplification arises for the laminar flows of in- 
terest here in that we can work in the low Mach number 
approximation: For the corresponding small flow veloc- 
ities the variations of the total pressure over the whole 
flow chamber domain will be much smaller than its ab- 
solute magnitude. Gradients of this small variation, de- 
noted as dynamic pressure p, can then be neglected in all 
equations except the momentum balance. 

Working in the low Mach number approximation the 
specific equations for the steady-state stagnation flow 
problem are e.g. discussed in detail by Kee, Coltrin and 
Glarborg 7 . The key assumption in the derivation is that 
variations of partial pressure, temperature and axial ve- 
locity in radial direction are much smaller than in axial 
direction, at least near the symmetry axis in the center 
of the catalyst. A lowest order expansion in the radial 
dependence leads then to a set of differential equations 
in which all fields depend only on the axial coordinate z. 



Here, A r = ^d r p = const, is the so-called radial pressure 
curvature, p the shear viscosity, k the thermal conduc- 
tivity, c p the specific heat capacity, c PjQ the specific heat 
capacity of species a, and j a ^ z the axial component of 
the diffusive mass flux of species a. 

Together with the ideal gas law, eq. (1), and in the 
low-Mach-number-limit the condition 



— —pkftT = const. 
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this set of equations allows to solve for all dependent 
fields, p(z), Y a (z), u(z), V(z), j a , z (z), T(z) andp(z). In 
fact, the axial momentum balance, eq. (4), is decoupled 
from the other equations. All fields except the dynamic 
pressure p can be determined without it, and it can there- 
fore be used afterwards to fix p from the other calculated 
fields. The problem gets fully determined by boundary 
conditions at the inlet and at the surface. Specifically, 
the second order equations demand independent infor- 
mation about V, T and Y a at both ends of the domain, 
and the first-order continuity equation requires informa- 
tion about u on one boundary. As this information about 
u is provided at both boundaries in the here discussed 
finite-gap stagnation flow, i.e. at inlet and surface, the 
resulting overdctcrmination is resolved by treating the 
unknown A r as an eigenvalue of the problem, i.e. its 
magnitude is adjusted to satisfy the additional boundary 
condition. 



2. Thermophysical and transport parameters 
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What remains for the numerical solution are values 
for the transport coefficients p and k, and a diffusion 
theory relating the diffusive mass flux j a „ to composition 
gradients. Further, we need expressions for the isobaric 
specific heat capacities c„ iQ and c p , as well as for the 
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TABLE I: Material parameters for the gas-phase species re- 
quired for the FD modeling: Characteristic diameter a and 
energy e, as well as vibrational frequencies uj and zero Kelvin 
component of the specific enthalpy h° a (which includes the 
zero-point energy contribution). The latter is referenced as 
usual with respect to the standard state of the atomic species, 
i.e. gas-phase oxygen and solid graphite for O and C, respec- 
tively. 
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CO 


C0 2 




a 


(A) 


3.458 


3.652 


3.769 


[8] 


e/k B 


(K) 


107.4 


98.1 


245.3 


[»} 


Hu) 


(meV) 


196 


269 


291, 167, 83 (2x) 


[9] 


h° 


(eV) 





-1.179 


-4.074 


[10] 



specific enthalpy h a (see section IIC below) . For the gas 
phase considered here the mixture specific heat is simply 
the mass-weighted sum of the species specific heats 



(8) 



For the c P;Q , themselves we consider translational, rota- 
tional and vibrational degrees of freedom, treating the 
prior two in the classical limit and the latter in the har- 
monic approximation 11 



with 
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, and N£ b and A^ ot the number of 
vibrational ana rotational degrees of freedom, respec- 
tively. For the three linear molecules O2, CO, and CO2 
jyrot _ 2 : anc } Table I lists their vibrational frequencies 
ui at i as taken from ref. 9. 

Computing the transport coefficients for a multi- 
component gas phase in complete generality is a very 
complex task of its own. For the accuracy level and 
gas-phase conditions of interest for our study effective 
semi-empirical molecular transport models provide for- 
tunately a fully sufficient description. From the manifold 
of such existing models we adopt a strategy described 
by Cloutman that is also used in the reactive flow CFD 
program COYOTE 8 . In a first step, this strategy relies 
on standard mixture-averaged approaches that relate fi 
and k of the multi-component gas mixture to the pure 
species values. In the case of the viscosity, this is the 
Wilkc formula 8 ^ 11 
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with 
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where m a is the mass of species a and fi a is the pure 
species viscosity. X a is the molar fraction of species a. 
For the thermal conductivity the analogue relation is 8 ' 11 
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(12) 



where the pure species conductivities. Note that 

$ Q( 3 is exactly the same as in Wilkes formula, i.e. it 
depends on the viscosities and not on the thermal con- 
ductivities. From the kinetic theory of dilute gases the 
expression 8 
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160F alnW*{T*) ' 
with T* = k B T/e a and 

(2,2)* (T *) = 1.147(7^-0.145 + (T * + 5 )2 t (w) 

provides the species viscosities in terms of two empirical 
parameters, a characteristic diameter a a and a character- 
istic energy e a . For a gas with Lennard- Jones interaction 
between the molecules a a and e a are the two parameters 
defining the interaction potential 7 . For the general case 
referencing to a Lennard- Jones model just provides a con- 
venient way of representing the temperature dependence 
of the transport coefficients and the two parameters need 
not to have a microscopic meaning. Values for these pa- 
rameters for a wide range of species are found in data 
bases and we list in Table I those for the species O2, CO, 
CO2 needed here 8 . From the thus determined viscosities 
the thermal conductivities are derived via the Eucken 
correction 8,11 



K (cVj.a: )Ma 
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(15) 



In the present application the remaining diffusive mass 
fluxes are predominantly driven by the concentration gra- 
dients and can then be implicitly obtained from corre- 
sponding Stefan-Maxwell equations 7 ' 11 at each point in 
the flow field 
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where the D°^S are the binary diffusion coefficients be- 
tween species a and /3. For their determination we em- 
ploy again an expression from the kinetic theory of dilute 
gases 8 ' 11 



n bin _ 



1/2 



rp3 ( m Q +mg \ 



1/2 



(17) 



2.1705 



with a a p = (cr Q + <Tp)/2, T* p = k B T / ^/e^ and 

tt (1A) *(T*p) = 1.0548(T* (3 )- ai5504 +(T^+0.55909)- 

(18) 

which also needs only the empirical characteristic diam- 
eter o a and energy e a of each species, cf. Table I. 
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C. Interface II: Surface 

As already mentioned, we need to specify a number of 
boundary conditions at the solid surface to fully deter- 
mine the set of stagnation flow equations. For the radial 
fluid velocity this is the standard no-slip condition, i.e. 
v(z = L) = 0. The normal component of the velocity at 
the surface (u(z = L) in our case) is commonly termed 
"Stefan velocity" 7 ' 12 . It is calculated by considering the 
mass balance at the surface 13 , and a finite Stefan veloc- 
ity accounts for transient storage or release of species due 
to a changing average surface composition. For the here 
discussed steady-state operation this is not the case and 
we have u(z = L) = 0. 

As the CO oxidation reactions at the surface are the 
only processes that consume reactants and yield products 
in stationary operation, the boundary condition for the 
diffusive mass flux is 

j*A z = L ) = -< url = -m a v a TO¥ , (19) 

where the chemical source term for each species is 
simply given by the overall rate of reaction events (per 
area and time) times the stoichiomctry coefficient v a and 
mass m a of the species in the reaction. For the simple 
CO oxidation reaction, vqo = — 1, vq, 2 = —1/2, and 
vco 2 = 1- 

The heat release connected to the exothermic conver- 
sions must also enter the heat balance at the surface. 
With the previously introduced boundary conditions the 
surface energy balance reduces to the mere requirement 
that the normal heat flux is continuous at the surface 

dT 

- k—(z = L) + ^h a j a ^ z (z = L) = q so iid,z(z = L) . 

a 

(20) 

Here, the two terms on the left hand side account for 
the heat transported away by the gas phase and the 
heat released by the surface chemical reactions, respec- 
tively, which must balance the heat flux into the solid 
<3soiid, z[z = L). The temperature-dependent contribu- 
tion to the specific enthalpies h a is determined com- 
pletely equivalcntly to the specific heat capacities, i.e. 
by considering translational, rotational and vibrational 
degrees of freedom as described above. The additionally 
required zero Kelvin component h° a including the vibra- 
tional zero point energies can be drawn from thermo- 
chemical tables 10 , and for completeness they are included 
in Table I for the species O2, CO, and CO2. 

At first glance this use of experimental thermochem- 
ical data might seem inconsistent with the use of first- 
principles energetics in the lp-kMC simulations. Our 
choice is motivated by the consideration that all other 
transport parameters are equally derived from experi- 
ment. The empirical transport models provide a fully 
sufficient and controlled description of the macro-scale 
flow structures at the accuracy level of interest for this 
study. As argued in more detail below such a description 
is reached for the reactive surface chemistry through the 



use of first-principles based lp-kMC modeling. We are 
thus faced with the problem of matching these two de- 
scriptions. In our approach this match occurs uniquely 
through the TOFs. Whatever energetics is required for 
their determination comes exclusively and consistently 
from first-principles, here density-functional theory with 
a semi-local exchange-correlation functional. Vice versa, 
the entire transport description is exclusively and con- 
sistently based on experimental numbers, such that the 
balance leading to its effective accuracy is not disturbed 
by occasional parameters as e.g. the h° a coming from ap- 
proximate first-principles theory. While we feel that such 
considerations are necessary for the envisioned error- 
controlled multiscale modeling, we note that in practice 
the semi-local functional employed in the lp-kMC over- 
estimates the zero Kelvin enthalpy change for the CO ox- 
idation reaction in Table I by only about 10%, such that 
none of the conclusions reported below would be touched 
when using this energetics instead of the experimental 
one. 

Since we are considering heat transport in both the 
solid and the gas phase, wc finally need two conditions de- 
scribing the change of the temperature field when cross- 
ing the surface. Equation (20) can only serve as boundary 
condition for the gas-phase heat transport. Another one 
is needed for the heat transport in the solid. Here, we 
want to assume that the temperature is continuous across 
the gas-solid interface. Within our interest in modeling 
the reactivity at near-ambient conditions, this should be 
rather well ensured by the frequent gas-surface collisions. 

D. Region II and Interface III: Sample and sample 
backside 

Neglecting possible heat-induced deformations or 
phase transformations we only account for heat conduc- 
tion through the single crystal. For the stagnation flow 
problem there is no radial variation of the gas-phase tem- 
perature. We maintain this description within the sample 
and correspondingly also model the heat transport as a 
one-dimensional problem described by Fourier's law 

dT 

<?solid,z = — Ksolid-T" ' (^1) 

dz 

where K SO Hd is the heat conductivity of the sample, which 
for simplicity we assume to be temperature-independent. 
Solution of this equation requires fixing a boundary value 
at the backside of the single crystal q S ohd,z(z = L + <i), 
which thus describes the degree of heat dissipation that 
is possible e.g. through radiative loss or the contact of 
the single crystal with the sample holder. Specifying this 
for real experimental setups is a demanding task and we 
suspect that in present in situ experiments this value will 
vary largely. Addressing these specificities in a quantita- 
tive way is clearly outside the capabilities of the present 
idealized reactor model. With real experimental setups 
lying anywhere in between we therefore analyze the rele- 
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vance of this factor for thin single crystals by considering 
two opposite extremes: First a fixed temperature at the 
sample backside to mimic a highly efficient heat coupling 
of the crystal to the system, and second a zero heat flux 
at the sample backside to represent a well insulated sam- 
ple. For either boundary condition eq. (21) can be solved 
analytically and provides in turn the missing boundary 
condition for the surface heat balance, eq. (20). 

Isolated sample (adiabatic limit) 

qsoUdA z = L ) = ■ ( 22 ) 
Connected sample (isothermal limit) 

fcoUd,.(* = £) = (T(z = L) - T back ) , (23) 

where T back is the temperature at the back of the crys- 
tal. While in principle this temperature could have any 
value (e.g. through controlled heating) we will assume 
here that it is identical to the one of the outside sys- 
tem and thus to the temperature of the incoming gas at 
the inlet, i.e. T back = T lnl . The remaining material pa- 
rameter to specify in the resulting boundary condition is 
the bulk heat conductivity. An experimental value for 
bulk Ru0 2 is K so iid(Ru0 2 )= 0.50 W cm" 1 K" 1 . 14 How- 
ever, almost all in situ work on RuO2(110) has in fact 
been performed on ultra-thin films grown on Ru crys- 
tals, which would suggest that the value K S oiid(Ru) = 1.17 
W cur 1 K _1 is more appropriate 15 . Fortunately, for the 
results reported below it makes no difference which value 
for this quantity is used. At the rather high thermal 
conductivity of either metallic Ru or the metallic oxide 
Ru0 2 this dissipation channel is so dominant, that - as 
soon as it is enabled in the surface heat balance, eq. (20) 
- it simply ensures that the surface temperature remains 
at the nominal value T surf (z = L) = T inl . Regardless of 
the specific value of K solld the second boundary condition 
is therefore simply equivalent to modeling an isothermal 
reactor limit, while the first boundary condition would 
correspond to the adiabatic limit. 

E. Numerical solution 

The stagnation flow equations, eqs. (1) and (2-7), can 
be transformed into a semi-explicit system of differential- 
algebraic equations (DAE) 

-rVi = f(Vj,C,z) , 

= g( yj ,C,z) , (24) 

where yj are the so-called differential components and 
C the algebraic component. In our case the differential 
components are density p, mass fractions Y a , temper- 
ature T, velocity components u and V, diffusive mass 



fluxes j atZ , intrinsic heat flux — K^-, pressure curvature 
A r , and The algebraic component is the gradient of 

the density, which is determined by the requirement that 
the total pressure is constant between inlet and surface. 

We solve the DAE boundary value problem numeri- 
cally using the COLDAE package 16 . This package uses 
piecewise orthogonal collocation at Gaussian points and 
has an adaptive mesh strategy allowing for an error- 
controlled solution. We use the default option controlling 
the number of intermediate points and an initial equidis- 
tant mesh with 10 spacings. In all simulations we employ 
a tolerance of 10 -4 for each differential component. In or- 
der to improve the stability and have full error control, all 
variables, dependent or independent, are presented in ap- 
propriate units. The employed units are the inlet-surface 
distance L for length, and for the velocity, density, and 
mass fractions their values at the inlet. The employed 
temperature scale is Kelvins. We use the representation 
(T — T lnl ) for the temperature, so that this renormalized 
temperature is always zero at the inlet. The mass fluxes, 
heat flux, and density gradient are expressed in multiples 

j-jcff.inlyinl in i inl ™ . . 

of L ° , ^jr-, and respectively, where D ctt - lnl is 

the mixture averaged diffusion coefficient 7 . Finally, the 
radial pressure curvature A r and p,?¥- are scaled with 

100p inl ^- and p ini ^, respectively. 

The software uses a (damped) Newton strategy to find 
a solution starting from an initial guess. The central 
features of the initial guess for the first simulation were 
constant p, T and Y a , as well third order polyno- 

mial for u that fulfills the boundary conditions. All re- 
maining unknowns were approximated according to these 
assumptions. For subsequent simulations we used the re- 
sults of previous simulations where appropriate. In those 
cases the adapted mesh from the previous simulation was 
coarsened to contain only half as many grid points as ini- 
tially. 



III. SURFACE REACTION CHEMISTRY: 
FIRST-PRINCIPLES KINETIC MONTE CARLO 
SIMULATIONS 

The actual surface catalytic activity enters the FD sim- 
ulations through the TOFs in the boundary value eq. 
(19) for the partial mass fluxes at the surface. The corre- 
sponding calculation of the rate of product formation per 
surface area and time from information about the elemen- 
tary processes in the catalytic cycle is the realm of mi- 
crokinctic theories. The most prominent such approach 
relies on phcnomenological rate equations which only 
consider the mean-field averaged concentrations (cover- 
ages) of the reaction intermediates at the surface 17 . This 
level of modeling of the surface chemistry is the prevalent 
standard in reactor engineering 7,12 . There, the kinetic 
quantities entering the rate expressions are in fact often 
treated as adjustable parameters. In a top-down fashion 
the idea is thus to use macroscopic reactor data to derive 
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effective insight into the on-going surface catalytic activ- 
ity. In more bottom-up oriented work the kinetic quan- 
tities are alternatively drawn from independent detailed 
experiments, or in modern hybrid approaches increas- 
ingly from first-principles calculations 18 ' 19 . The idea of 
such integrated approaches is correspondingly to model 
how both the intrinsic surface chemistry and transport 
effects contribute together to the macroscopically observ- 
able activity in a given reactor setup. 

The latter is also the central objective of the present 
study. However, for the here aspired quantitative mod- 
eling present-day hybrid approaches are not sufficient. 
Use of scattered experimental and first-principles kinetic 
data from different sources and in the latter case poten- 
tially from different levels of approximate theory incurs 
a rather uncontrollable error. Even if the kinetic pa- 
rameters of all involved elementary processes were reli- 
able, there is still the error from the approximate mean- 
field treatment underlying the rate equation approach. 
In fact, for the here studied CO oxidation reaction at 
RuO2(110) this error has recently been shown to be qual- 
itative with deviations of the mean-field TOFs spanning 
up to several orders of magnitude 20 . Aiming at an error- 
controlled multiscale modeling of predictive quality we 
therefore employ for the description of the surface ki- 
netics lp-kMC as most accurate approach with explicit 
account of the correlations, fluctuations and detailed spa- 
tial distributions of the chemicals at the surface. 



A. lp-kMC model of CO oxidation at RuO 2 (110) 

The molecular-level basis of lp-kMC is a microscop- 
ically correct first-principles description of the elemen- 
tary processes involved in the catalytic cycle. In the es- 
tablished model of CO oxidation at RuO 2 (H0) 5 ' 6 this 
is specifically the set of 26 elementary processes defined 
by all non-correlated site and element specific adsorp- 
tion, desorption, diffusion and reaction events that can 
occur on a lattice spanned by two different active sites 
offered by the surface, so-called bridge (br) and coordi- 
natively unsaturated (cus) sites. For all these processes 
density-functional theory in conjunction with harmonic 
transition-state theory is used to compute the kinetic pa- 
rameters. The resulting 26 first-principles rate constants 
form then the essential input for the actual lp-kMC sim- 
ulations which evaluate the statistical interplay among 
the surface chemical processes by following the long-term 
time evolution of the open catalytic system through nu- 
merical solution of a Markovian master equation 1 . 

Using exactly the computational setup as detailed 
before 5 ' 6,20 , these simulations are carried out for the 
present purposes for a given local temperature T surf and 
reactant partial pressures Po 2 rf and p^Q directly at the 
surface, which in particular fixes the impingement and 
therewith the rate constants of the adsorption processes. 
Under such conditions, the system eventually reaches a 
unique steady-state, in which the detailed surface compo- 



sition and occurrence of the individual elementary pro- 
cesses still exhibit the correct microscopic fluctuations, 
yet when averaged over simulation cells exceeding the 
characteristic correlation lengths at the surface have well- 
defined and constant values. These values thus comprise 
the average rate of reactant to product conversion under 
the given gas-phase impingement and local temperature, 
i.e. exactly the TOFs that enter the partial mass flux 
boundary condition of cq. (19). While it is only this 
averaged quantity that matters for the macroscopically 
described flow field, it is important to note that it is still 
properly derived from microscopic simulations that fully 
account for the site heterogeneity and distributions at 
the surface. This is thus distinctly different to the men- 
tioned mean-field based phenomenological descriptions 
that are commonly integrated in the CFD modeling of 
macro-scale flow structures. 



B. Integration of lp-kMC into the fluid dynamical 
environment 

The lp-kMC and FD simulations are intricately cou- 
pled. On the one hand, the TOFs required in eq. (19) 
to close the stagnation flow equations are provided by 
the lp-kMC simulations. On the other hand, fixing the 
surface impingement in the lp-kMC simulations requires 
the local temperature and gas-phase partial pressures di- 
rectly at the surface, which are determined by the heat 
and mass transport modeled at the continuum level. A 
straightforward approach to this interdependence is a si- 
multaneous solution until self-consistency between flow 
and lp-kMC boundary condition is achieved 21 . For the 
here discussed stagnation flow equations this approach 
is in principle feasible 22 , albeit potentially numerically 
unstable 23 . However, for more complex reactor geome- 
tries it would quickly become intractable, as usually sev- 
eral independent lp-kMC simulations would be required 
for every spatially resolved cell at the surface. Precisely 
due to the necessity to continuously rerun lp-kMC simu- 
lations the approach would also be very inefficient for the 
here intended simulation of catalytic activity at a large 
variety of flow conditions. 

We instead achieve a computationally much more ef- 
ficient formulation by decoupling the interdependence 
through an instantaneous steady-state approximation 12 . 
The kMC simulations are first carried out to determine 
the steady-state TOFs for a wide range of surface im- 
pingement and local temperature conditions. The result- 
ing grid data is then interpolated to a continuous rep- 
resentation, which in turn provides the entire necessary 
boundary condition for the stagnation flow problem. For 
the steady-state operation targeted in this manuscript 
this divide-and-conquer type approach is exact and may 
easily be applied to more complex reactor geometries. It 
could even be extended to transient phenomena under 
the assumption that on the time scale of relevant flux 
variations the surface chemistry adapts quasi instanta- 
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neously to the new steady-state, hence the name. 

For the CO oxidation at RuO 2 (110) we thus first com- 
puted lp-kMC steady-state TOFs for the entire relevant 
range of temperatures and gas-phase composition. With 
a negligible CO2 readsorption probability these TOFs are 
independent of the CO2 partial pressure. Specifically we 
then used a dense grid with 25 K spacing for the tem- 
perature range 400 K < T surf < 850 K and with logscale 
spacing to cover the range 10 -6 atm < Po" f < 10 2 a ^ m 
with 30 and the range 10" 5 atm < p s gg < 10 2 atm 
with 42 spacings. Through modified quadratic Shepard 
interpolation 24 ' 25 this is converted into a reliable contin- 
uous representation TOF(T surf ,pg l 2 rf ,p s c u <5 f ) that is finally 
presented as boundary condition to the stagnation flow 
solver. 



IV. RESULTS 

The intrinsic activity resulting from the lp-kMC model 
for the CO oxidation at RuO2(110) has been analyzed 
for a wide range of temperature and pressure conditions 
before 5,6 . Not surprisingly, high catalytic activity is only 
observed for a rather narrow range of gas-phase condi- 
tions, which stabilize O and CO simultaneously at the 
surface in appreciable amounts. For more O-rich feed 
the surface is poisoned by oxygen, for more CO-rich feed 
the surface is poisoned by CO, and little CO2 is formed in 
cither case. For gas-phase conditions in the UHV regime, 
which allow direct measurements of the intrinsic activity, 
the model reproduces existing experimental TOF data 
quantitatively 5,6 . It must be stressed though that for 
the description of the CO-poisoned regime the model has 
clear limitations. In corresponding CO-rich feeds the ox- 
ide surface would in reality eventually be reduced. By 
construction, this and the catalytic activity connected 
to such a phase transformation can not be grasped by 
the present lp-kMC model assuming an intact underly- 
ing RuO 2 (110) lattice. 

Notwithstanding, this restriction concerns the model- 
ing of the surface chemistry. The focus of the present 
work is instead to quantitatively integrate a given mi- 
crokinctic description of this surface chemistry into a FD 
framework to assess heat and mass transport effects in 
a reactor geometry representative for in situ studies of 
model catalysts. For this endeavor the existing account 
of the surface chemistry is - despite its noted limitation 
- very suitable, in particular as it exhibits a number of 
features that we consider rather generic for a high-TOF 
reaction like the CO oxidation: (i) The intrinsic catalytic 
activity is narrowly peaked in a small range of gas-phase 
conditions, (ii) This activity is not sufficiently described 
by standard rate equation formulations. For predictive 
quality the first-principles based microkinetic modeling 
must therefore be based on an approach like lp-kMC 
that explicitly accounts for the detailed spatial distribu- 
tion of the chemicals at the surface. In turn, it is the 
latter type of approach to which the FD environment 



must be coupled, e.g. through the instantaneous steady- 
state approximation employed here, (iii) The peak activ- 
ity at optimum partial pressures increases rapidly in the 
temperature range of interest. Towards the upper end at 
around 500-600 K and together with the high exothermic- 
ity of the CO oxidation reaction, this leads to a degree 
of mass conversion and heat release at the surface that is 
prone to transport limitations in the reactor. 

As we will see the amount of heat dissipation possible 
at the back of the thin single-crystal is a crucial factor 
for such limitations that can mask the true intrinsic ac- 
tivity in in situ model catalyst studies. Not aiming (nor 
being able to) give a detailed account for one specific ex- 
perimental setup we will analyze this in the following in 
more generic terms for the two already described oppo- 
site limits. In the adiabatic limit there is no heat flux 
at all through the sample backside, mimicking to some 
degree the situation that could e.g. arise from an insu- 
lating sample holder. In contrast, in the isothermal limit 
we assume the sample to be sufficiently well connected to 
the outside system that a constant nominal temperature 
is maintained throughout. Here, this is chosen to be the 
same temperature as that of the gases at the inlet. 

Considering the peaked structure of the intrinsic cat- 
alytic activity in (T,Pq" 1 ,p^q) -space we can conve- 
niently study heat and mass transfer effects in these two 
limits focusing on prototypical sets of gas-phase con- 
ditions: For defined temperature, essentially zero CO2 
concentration (j>co = 10~ 5 atm throughout) and near- 
ambient oxygen partial pressure at the inlet these sets 
comprise a range of inlet CO partial pressures. They 
cover the O-poisoned regime at the lowest Pqq , the CO- 
poisoned regime at the highest Pqq, and span over the 
conditions of highest intrinsic activity. 



A. Adiabatic limit 

1. Surface heating 

The heat flux through the back of the sample is com- 
pletely suppressed in the adiabatic limit. The only dissi- 
pation channel left for the heat released by the exother- 
mic surface reactions is then into the surrounding gas 
phase. Compared to the heat conduction through a 
metallic sample this channel is rather inefficient. One can 
therefore suspect that it may not be efficient enough to 
maintain the nominal surface temperature, once the TOF 
and therewith the generated heat rate exceeds a certain 
critical value. Instead, the surface will heat up and give 
rise to gas-phase temperature gradients from inlet to ac- 
tive surface. This is indeed what we find in the coupled 
lp-kMC-FD simulations. When using representative pa- 
rameters for the inlet distance L = 1 cm and axial inlet 
flow velocity u lnl = lem/s significant deviations from 
the nominal surface temperature set in for near-ambient 
environments at TOFs exceeding ~ lOsite -1 s" 1 . Such 
peak TOFs are reached for optimum partial pressure con- 
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FIG. 2: (Color online) Comparison of intrinsic steady-state 
TOFs as resulting from lp-kMC (black thin dashed line) 
with observable TOFs when accounting for transport effects 
in the stagnation flow reactor in the adiabatic limit. For 
T lnl = 500 K and pq\ = 0.3 atm the shown range of inlet 
CO partial pressures spans from the O-poisoned to the CO- 
poisoned regime with the intrinsic "most active" state reached 
at intermediate pco corresponding roughly to a stoichiomet- 
ric feed. The suppressed heat flux at the sample backside 
allows the system to sustain a high-activity operation mode 
for more CO-rich conditions than this nominal most active 
state. In this high-activity branch of the TOF-profile the 
surface temperature, J nsurf shown in the lower panel, is signif- 
icantly increased. The activity and extension of the branch 
depends on the reactor details. Shown here are results for 
constant inlet velocity u ml = lcm/s and varying inlet-surface 
distance: L — 1mm (dotted red line), L = 1cm (solid red 
line) and L = 10 cm (dashed red line). 



ditions at temperatures above about 500 K. 

The effect of the ensuing temperature increase at the 
active surface on the observable steady-state conversion 
rates is quite dramatic already at this threshold tempera- 
ture and is illustrated in Fig. 2. While the intrinsic activ- 
ity is peaked close to stoichiometric feeds, the heat trans- 
port limitations lead to the stabilization of a high-activity 
operation mode that extends to significantly more CO- 
rich conditions. In this branch of the observable TOF- 
profile shown in Fig. 2 the surface temperature is up 
to 150 K higher than the nominal temperature at the in- 
let. Quantitatively this and the concomitant TOFs in 



the newly established high-activity mode depend on the 
details of the reactor setup. This is exemplified in Fig. 
2 with the TOF-profiles that result when increasing or 
decreasing the inlet distance by one order of magnitude. 
For a larger L = 10 cm the extension of the high-activity 
branch is reduced at overall lower conversion rates, while 
for a smaller L = 1 mm the branch extends to much 
higher CO partial pressures and the observable conver- 
sion rate exceeds the peak intrinsic activity by more 
than one order of magnitude. Similar, but quantitatively 
smaller variations are obtained when changing the in- 
let velocity by one order of magnitude up or down (not 
shown). For a lower w ml = lmm/s the high- activity 
branch exhibits slightly smaller TOFs and extends up to 
slightly smaller p^Q than for the w ml = lcm/s displayed 
in Fig. 2. An increase to u = lOcm/s, on the other 
hand, increases the extension of the branch to higher CO 
partial pressures at also higher TOFs than those shown 
in Fig. 2. 

Regardless of these quantitative variations, the net ef- 
fect of the surface heating resulting from the transport 
limitations is in all cases a substantially changed observ- 
able TOF-profile compared to the true underlying intrin- 
sic reactivity. The absolute TOFs in the relevant high- 
activity regime are significantly different and the inlet 
gas-phase conditions for which highest conversions are 
obtained are shifted to much more CO-rich feeds. Ob- 
viously, if these observable TOFs were mistaken for the 
intrinsic reactivity, wrong mechanistic conclusions about 
the on-going surface chemistry in such in situ environ- 
ments would be derived. Furthermore, the observable 
TOF-profile in Fig. 2 exhibits another feature that is 
completely absent in the intrinsic reactivity: For a range 
of CO partial pressures we obtain two stationary solu- 
tions: The high-activity branch and in addition a low- 
activity branch that coincides with the intrinsic activity. 
As the underlying lp-kMC model of CO oxidation at 
RuO2(H0) has no multiple steady-states 20 , this bistabil- 
ity arises solely from the coupling of macroscopic trans- 
port and surface chemistry. 



2. Formation of a boundary layer 

An analysis of the observed range of transport effects 
(high-activity branch with concomitant bistability and 
variations with reactor setup) starts from the anticipated 
formation of a finite boundary layer above the flat-faced 
model catalyst. As schematically drawn in Fig. 3 non- 
vanishing gradients of temperature and partial pressures 
are in general restricted to this boundary layer. This is 
due to the convective nature of the transport. The gas 
streams towards the surface and the chemical reactions at 
the surface have no influence on the fields far away. Heat 
and species move with the flow, which dominates over any 
non-convective transport. Only when the axial velocity 
approaches zero near the surface, diffusion and heat con- 
duction kick in, as they are now of similar size as the 
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FIG. 3: (Color online) Illustration of the temperature and 
pressure profiles in steady-state stagnation flow. Diffusion 
and heat conduction only take place in a boundary layer above 
the surface; outside transport is purely convective. In the here 
discussed adiabatic limit sufficiently large TOFs lead to an 
increase of temperature in this boundary layer, while in the 
isothermal limit they lead to a change of partial pressures. 
In the general case both effects can be intricately intermin- 
gled. The boundary layer expands when the inlet velocity is 
decreased, ultimately filling the whole gap between inlet and 
surface in the limit n lnl — > 0. The variation with the inlet- 
surface distance L is opposite, i.e. the boundary layer shrinks 
with decreasing L. 

convective transport. For the here discussed adiabatic 
limit the continuous catalytic conversions at the surface 
lead to a continuous heat release into the gas phase. For 
sufficiently large TOFs above the critical value this heat 
rate is too high to be efficiently transported away by the 
flow. The result is a steady-state temperature profile 
above the surface as sketched in Fig. 3, where the tem- 
perature rises within the boundary layer from T ml to the 
values for T su shown in Fig. 2. On the other hand, sur- 
face mass conversions of the order of magnitude as those 
of Fig. 2 arc still too low to significantly affect the gas 
composition in the boundary layer. The formed products 
are transported away sufficiently quickly so that the par- 
tial pressures remain essentially unchanged, i.e. we find 
psurf ^rni f or gjj conditions discussed in Fig. 2. 

The latter feature allows us to suitably discuss the ob- 
served transport effects in terms of the intrinsic reactivity 
summarized in Fig. 4. Shown is the TOF profile as ob- 
tained with lp-kMC for the same p l §\ w p s ^ f = 0.3 atm 
as in Fig. 2 and as a function of surface temperature and 
CO partial pressure. The region of highest-activity is 
narrowly peaked and the peak TOFs at this rim increase 
steadily with temperature. A crucial feature for the fol- 
lowing is that the rim does not go along constant partial 
pressure ratio for higher temperatures (i.e. vertically up 
in Fig. 4), but shifts continuously towards more CO-rich 
mixtures (i.e. diagonally up in Fig. 4). This is because 
the most active state is characterized by the coexistence 
of O and CO at the surface, which follows more a con- 
stant ratio of chemical potentials than partial pressures. 
Imagine now that the system starts at the nominal inlet 
and surface temperature of 500 K, corresponding to the 
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FIG. 4: (Color online) Intrinsic TOF contour plot as com- 
puted with lp-kMC for the same constant p 1 ^ ~ pSj = 
0.3 atm as in Fig. 2. Marked by red lines are the TOF and 
surface temperature that result in the steady-state stagnation 
flow with it = lem/s and varying inlet distances, i.e. the 
conditions behind the different red lines in Fig. 2. The line 
shapes follow the ones of Fig. 2, i.e. L = 1 mm (dotted red 
line), L = 1cm (solid red line) and L = 10cm (dashed red 
line). 

bottom horizontal line in Fig. 4. As long as the intrin- 
sic TOF does not exceed the critical value, the system 
is able to maintain this temperature and the observable 
TOF is identical to the intrinsic one, cf. Fig. 2. For the 
range 0.4 atm < Pqq < 0.85 atm the activity is, however, 
above the critical value. The system cannot transport 
the generated heat away sufficiently quickly, and surface 
and boundary layer start to heat up. As directly appar- 
ent by focusing on e.g. p™Q = 0.5 atm in this critical 
pressure range the intrinsic TOFs increase with increas- 
ing surface temperature (vertically up in Fig. 4). This 
gives rise to a runaway process. Higher TOFs generate 
more heat, which increases the surface temperature and 
therewith leads to even higher TOFs. This will only stop 
once the system is over the highest activity rim in Fig. 
4. Further increase in surface temperature (i.e. moving 
even more vertically up along the line of constant Pqq 
in Fig. 4) leads then to a decrease in intrinsic reactivity. 
This allows the system to find a new steady-state, viz. 
the high-activity branch corresponding to the increased 
surface temperature in Fig. 2. 

Where on this downward TOF slope the system pre- 
cisely settles down depends on how efficiently the gen- 
erated heat can be transported away, i.e. which sur- 
face temperature results for the heat rate connected with 
a given intrinsic TOF. This is crucially controlled by 
the thickness of the boundary layer and corresponding 
thickness variations rationalize the entire observed de- 
pendence on the reactor setup. A smaller inlet distance 
compresses the boundary layer and therewith enables a 
better heat convection. Correspondingly, at smaller L 
the system ends up closer to the most active rim in Fig. 4 
and the high-activity operation mode exhibits higher ob- 
servable TOFs as seen in Fig. 2. A higher inlet velocity 



11 



has the same effect on the boundary layer, and therewith 
also leads to higher observable TOFs as described above. 



The intrinsic TOF profile shown in Fig. 4 also helps 
to understand why the high-activity branches eventually 
break down and how their extension varies with the re- 
actor setup. Apart from the highest-activity rim domi- 
nated by the oxidation reactions between CO and O ad- 
sorbed at the most active cus sites one can discern at 
the upper edge of Fig. 4 a weak new rise of the intrin- 
sic TOFs. At these highest temperatures shown this is 
due to the "entropic" widening of the kinetic "phase" 
transition region where the appropriate chemical poten- 
tial ratio in the gas-phase stabilizes the coexistence of 
both reactants at the surface 6,26 . When over the rim 
the hitherto described decrease of the intrinsic TOF with 
increasing temperature will therefore eventually change 
into a new increase. If the system reaches this change 
of slope, a new runaway cycle of increasing temperature 
and TOF will start and no stationary operation mode 
can be stabilized (at least not in the temperature range 
of interest for the present study). Correspondingly, in 
Fig. 4 the high-activity branches always break down at 
positions of such a gradient change in the intrinsic TOF 
profile. As apparent from Fig. 4 the further away from 
the rim the high-activity branch is situated, the earlier it 
hits this slope change, i.e. its extension reaches only up 
to smaller Pqq. In the present adiabatic limit, modifi- 
cations of the reactor setup that compress the boundary 
layer (either by higher axial inlet velocity or smaller in- 
let distance) lead therefore to higher absolute observable 
TOFs in a high-activity branch that extends up to higher 
CO partial pressures. 



Finally, some remarks about the observed bistability 
are at place. The structure of the TOF profile in Fig. 4 
rationalizes why for some pressure conditions two steady- 
state solutions are obtained. One, in which the system 
exhibits a low activity that coincides with the intrin- 
sic one, and one, in which significant surface heating 
has brought the system above the highest-activity rim. 
While intuitive, the rationalization in terms of thermal 
runaway is at present clearly an interpretation. A verifi- 
cation would require the extension of the present steady- 
state approach to transient operation. Only correspond- 
ing time-resolved simulations will then give access to the 
wealth of phenomena that are now only suggested by the 
observed bistability. Notably this is the possibility of 
oscillations between the two modes. In contrast to e.g. 
purely surface reaction-surface diffusion driven oscilla- 
tions on single-crystals in UHV 27 the mechanistic details 
behind corresponding reactor-reaction oscillations in the 
ambient pressure regime are only poorly understood 28 . 
Obviously, extending the present model in this direction 
offers the prospect of a detailed analysis, on which we 
will concentrate in future work. 
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FIG. 5: (Color online) Comparison of intrinsic steady-state 
TOFs as resulting from lp-kMC (black thin dashed line) with 
observable TOFs when accounting for transport effects in the 
stagnation flow reactor in the isothermal limit. As in Fig. 2 
the data is for p'oj = 0.3 atm, but now for a higher temper- 
ature T inl = 600 K. The higher intrinsic peak TOFs at this 
higher temperature are significantly masked by mass trans- 
fer limitations. The variation with the reactor setup is il- 
lustrated for constant inlet velocity u m = lem/s by varying 
inlet-surface distances, L = 1mm (dotted red line), L — 1cm 
(solid red line) and L = 10cm (dashed red line). Using the 
same line styles the lower panel shows how the partial pres- 
sure ratio at the surface deviates from the nominal one at the 
inlet. 



B. Isothermal limit 

In the opposite isothermal limit the high thermal con- 
ductivity of the metallic sample allows for such an effi- 
cient removal of the generated heat that even at higher 
temperatures around 600 K, where the intrinsic peak 
TOFs at optimum partial pressures exceed 10 4 site _1 s _1 , 
no significant surface heating results. The temperature 
remains at the nominal inlet value throughout the en- 
tire system. As shown in Fig. 5 the intrinsic activity is 
nevertheless noticeably masked, this time by mass trans- 
fer limitations in the boundary layer. At it lnl = lem/s 
and L = 1 cm the maximum observable TOFs are lower 
than the peak intrinsic ones, and a high-activity branch 
extends now to much more CO-poor feeds. As in the 
adiabatic limit there is a range of CO partial pressures 
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for which we observe a bistability, and the results depend 
again quantitatively on the reactor setup. This is illus- 
trated in Fig. 5 by showing the data also for an increased 
(L = 10 cm) and decreased (L = 1 mm) inlet distance. 
For smaller inlet distances higher observable TOFs result. 
Highly comparable variations are obtained when chang- 
ing the axial inlet velocity by one order of magnitude 
up or down, with smaller it lnl yielding larger observable 
TOFs (not shown). The varying reactor conditions also 
affect the extension of the high-activity branch, which 
at maximum reaches down to j»J?q = 0.6 atm, i.e with 
Pol = 0.3 atm to stoichiometric feed. In addition, there 
is now in principle also a dependence of the results on 
the employed sample thickness d, which enters through 
the boundary condition, eq. (23). However, due to the 
high thermal conductivity this dependence is in practice 
negligible. Compared to the results in Fig. 5, which were 
obtained using d = 1 mm, changing the thickness by one 
order of magnitude up or down has virtually no effect on 
the observable catalytic activity. 



1. Mass transfer limitations 

The entire range of transport effects observed in the 
isothermal limit is this time due to mass transfer limita- 
tions in the boundary layer. At the high intrinsic TOFs 
around peak performance the mass conversion at the ac- 
tive surface is so large that these limitations lead to no- 
ticeable changes of the partial pressure profiles from inlet 
to surface. As schematically shown in Fig. 3 there is es- 
sentially a build-up of a significant product concentration 
that is no longer sufficiently quickly removed. This goes 
hand in hand, cf. eq. (7), with a decrease in reactant 
partial pressures, i.e. O2 and CO are hindered in their 
access to the active surface. Due to the similar trans- 
port parameters and mass of both diatomic reactants, cf. 
Table I, this limitation affects both species similarly. As 
long as the nominal inlet composition is different from 
stoichiometric feed, a corresponding roughly equal re- 
duction of both reactant partial pressures close to the 
surface will then effectively change the Pqo /Pol ra ti°- 
As also apparent from Fig. 4 at T lnl = 600 K the range 
of peak intrinsic activity corresponds to quite CO-rich 
feeds. In this range the mass transfer limitations there- 
fore lead to a noticeable increase of the Pco iPof ratio 
compared to the nominal inlet composition as shown in 
Fig. 5. At a nominal inlet composition that would cor- 
respond to optimum intrinsic activity, Pqq w 3 atm in 
Fig. 5, the surface then sees a comparatively more CO- 
rich feed and the observable TOFs are lowered compared 
to the intrinsic ones. On the other hand, at a nomi- 
nal inlet composition only slightly more CO-rich than 
stoichiometric feed, where the intrinsic activity would al- 
ready have collapsed in Fig. 5, the significantly more 
CO-rich feed effectively seen by the surface corresponds 
in fact to conditions close to optimum intrinsic activity. 
The observable TOF is much increased and the high- 
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FIG. 6: (Color online) Same as Fig. 5, but including the upper 
TOF limit set entirely by mass transfer (dotted line), see text. 
The blue dotted line indicates the range limited by oxygen 
mass transfer, the green dotted line the range limited by CO 
mass transfer. Shown is data for Pq\ — 0.3 atm, T lnl = 600 K, 
u' nl = lem/s, L = 1cm. 



activity branch of Fig. 5 results. Exactly at stoichiomet- 
ric feed this effective CO enrichment close to the surface 
ends, and for even more CO-poor mixtures possible mass 
transfer limitations would rather suppress the CO mi- 
nority species. However, for such partial pressure ratios 
the intrinsic activity is low anyway and no mass transfer 
limitations arise. At the latest the high-activity branch 
therefore breaks down at stoichiometric feed. With this 
understanding the observed variations with the reactor 
setup are also easy to rationalize. A smaller boundary 
layer as resulting from increased axial inlet velocity or 
reduced inlet distance reduces the mass transfer limita- 
tions. The partial pressure ratio at the surface gets closer 
to the nominal one. In turn, the observable TOFs ap- 
proach the intrinsic ones as illustrated in Fig. 5 for the 
varying inlet distances. 

In the presence of such mass transfer limitations a nat- 
ural question is to what degree they mask the intrinsic 
surface reactivity. Is the observable TOF profile the re- 
sult of a complex mixture of the on-going surface chem- 
istry and the gas-phase transport, or are the flow con- 
ditions in the reactor such that the measured activity 
conveys little information about the actual catalyst any- 
more. To qualify this it is useful to assess the upper 
TOF limit that results if mass flux is the only limitation. 
Such an estimate can be obtained by realizing that the 
steady-state mass conversion by the catalyst can never be 
higher than the one that completely depletes the minor- 
ity species at the surface. Rather than using the catalyst 
specific boundary condition eq. (19) that depends on the 
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actual intrinsic TOFs, we then simply employ for O-rich 
feeds 

pg£ = , (25) 

and for C O-rich feeds 

Po 2 rf = ■ (26) 

For the respective majority species the boundary condi- 
tion is still cq. (19), but the TOF entering this equation is 
now determined by the mass flux for the minority species, 
i.e. the conversion is completely dictated by the amount 
of impinging minority species. Figure 6 shows the upper 
TOF limit that results from this estimate for the afore 
discussed gas-phase conditions of Fig. 5. Apparently, 
the observable TOFs come very close to this upper limit 
for most of the active region. This means that in this 
regime the measurable profile has very little to do with 
the actual RuO2(110) catalyst, it rather reflects only the 
flow conditions in the employed reactor. Obviously, and 
similar to the adiabatic limit discussed before, if corre- 
sponding effects are not appropriately accounted for in in 
situ studies, wrong conclusions about the surface chem- 
istry at technologically relevant gas-phase conditions will 
be derived. 



V. SUMMARY AND CONCLUSIONS 

In summary we have presented an efficient approach 
to couple first-principles kinetic Monte Carlo and fluid 
dynamical simulations in the context of heterogeneous 
catalysis. This augments the accurate description of the 
surface chemistry achieved by lp-kMC with a continuum 
account of the heat and mass transfer in a given reactor 
geometry, or vice versa it integrates the accurate lp-kMC 
microkinetics into reactor-level modeling. In prevalent 
chemical engineering approaches the latter modeling in- 
stead incorporates phenomcnological microkinetic treat- 
ments based on mean-field rate equations. In contrast to 
such a description, the presented lp-kMC based multi- 
scale modeling approach derives the average flux quanti- 
ties required for the macroscopically described flow field 
properly from microscopic simulations that fully account 
for the site heterogeneity and distributions at the cata- 
lyst surface. As such it has the potential to carry the 
predictive power of the underlying electronic structure 
calculations for the elementary processes all the way to 
the reactor level. 

On the way to such a first-principles chemical engineer- 
ing we have applied the approach to the problem of in situ 
studies of model catalysts using the ambient pressure CO 
oxidation at RuO2(110) as a representative example. As 
a suitable, though idealized reactor model to discuss the 
transport at the flat-faced surface we have chosen a stag- 
nation flow geometry The observed catalytic function 
depends sensitively on the employed reactor geometry 



(mimicked in this study by varying the inlet-surface dis- 
tance) and the applied throughput rate (i.e. the stream- 
ing velocity at the inlet). For the thin, well heat conduct- 
ing single crystal the degree of heat dissipation possible 
at the back of the sample (e.g. through radiative loss or 
contact to the sample holder) is a further crucial factor. 
Not aiming, nor being able to address specific experimen- 
tal realizations this was considered through two opposite 
extremes: The isothermal limit mimicking a highly ef- 
ficient heat coupling of the crystal to the system, and 
the adiabatic limit to represent a well insulated sam- 
ple. In both limits transport effects were found to readily 
mask the intrinsic catalytic function at the high conver- 
sion rates reached at near-ambient gas-phase conditions. 
In the adiabatic limit this is due to a significant sur- 
face heating, in the isothermal limit due to mass trans- 
fer limitations that lead to the build-up of a significant 
concentration of products in the boundary layer above 
the active surface. With the single-crystal in real experi- 
mental setups neither perfectly heat-coupled nor isolated, 
these two effects discussed here separately will obviously 
be intricately intermingled and need to be disentangled 
by dedicated measurements and setups. Furthermore, 
we obtained in both limits a range of gas-phase condi- 
tions where the system exhibits two stationary operation 
modes, a low-activity branch corresponding to the in- 
trinsic reactivity and a high-activity branch which arises 
from the coupling of the surface chemistry to the sur- 
rounding flow field. A corresponding bistability obtained 
here in the steady-state limit clearly suggests that the 
system could oscillate between the two modes, possibly 
even inhomogeneously in form of reaction fronts over the 
single-crystal surface. In case of heat transfer limita- 
tions, an intuitive propagation mechanism would hereby 
be via the formation of local hot spots, while in the mass 
transfer case it would be via gas-phase coupling, with the 
presented approach establishing the intriguing possibility 
to quantify these model conceptions with first-principles 
based simulations. 

The main objective of in situ studies of model cata- 
lysts is a detailed, atomic-scale analysis of the catalytic 
function at technologically relevant gas-phase conditions, 
thereby bridging the pressure gap to the at present often 
much better characterized function in UHV. The range 
of transport effects discussed in the present study quali- 
fies the additional complexity that needs to be accounted 
for in corresponding work to prevent wrong mechanistic 
conclusions about the surface chemistry at high pressure. 
That this complexity has potentially not yet been suffi- 
ciently appreciated may very well be the reason for the 
many existing controversies in the field. Also because of 
the limitation of the employed lp-kMC RuO2(110) model 
with respect to a reduction of the oxide catalyst we have 
refrained from comparing our simulations to already pub- 
lished experimental data. Nevertheless, we note that the 
gas-phase conditions and TOFs discussed here are of the 
order of magnitude presented in a manifold of in situ 
studies of CO oxidation at late transition metal cata- 
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lysts. In this respect it is important to recognize that the 
CO oxidation reaction - that has been a fruitfly reaction 
in UHV Surface Science due to its alleged "simplicity" 
and model character - requires particular attention. The 
high turnovers that can be reached precisely because of 
this " simplicity" make this reaction much more prone to 
transport effects than other more complex, selective ones. 
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